从密度矩阵产生自然轨道-理论篇
1.自然轨道的定义
对于一个单或多行列式波函数方法（例如RHF, MP2, CCSD, CASCI, CASSCF等等），可将电荷密度（charge density）展开到一组正交归一的轨道上
注意电荷密度是个全空间的函数。其中nmo表示分子轨道数，是个厄米矩阵（实数下即实对称阵），此处表示分子轨道基下的密度矩阵，矩阵的迹必须等于体系电子数
这里所说的满足正交归一关系是指
其中为原子基函数，nbf为基函数数目，为原子基函数之间的重叠积分。
对于一个给定的体系和确定的波函数方法，是固定的，因此若换成另一组正交归一的轨道 ，便会对应一个新的矩阵，写成公式就是
其中
这其中有个特殊的酉变换尤其重要：存在一个特殊的使得是对角矩阵，即
此即为自然轨道表象，对应的分子轨道称为自然轨道（natural orbital, NO），不妨专门记为。此时可简写为
即为，称为第个自然轨道的占据数（natural orbital occupation number, NOON）。这式子很好理解：用电子数作为权重乘以轨道模方。所有轨道占据数加起来即为体系总电子数
举个简单的例子，在RHF方法里取值只能是整数2/0，对应双占据/空轨道
电荷密度简化为
即求和指标从所有轨道（nmo）减小为双占据轨道（occ），接着将分子轨道展开至原子基函数（AO basis）上，便可出现大家在量子化学课上熟知的RHF（原子基）密度矩阵元，写成矩阵形式
如果有一组轨道是当前轨道的酉变换，密度矩阵不会变，占据轨道的占据数也仍是2，即
因此RHF正则占据轨道（及其任意酉变换）本身也是自然轨道。
对于UHF则有四种常见的密度矩阵：alpha自旋，beta自旋，自旋密度矩阵（即alpha-beta密度矩阵差），总密度（即alpha+beta密度矩阵和），对应四种自然轨道：alpha自然轨道，beta自然轨道，自旋自然轨道(SNO)，UHF自然轨道（UNO）。前两种的轨道占据数严格为1，因此alpha/beta正则占据轨道（及其任意酉变换）本身亦是自然轨道；后两种则需要对角化相应的密度矩阵得到（见下文），轨道占据数是小数，SNO占据数范围[-1,1]，UNO占据数范围[0,2]。
在一般的方法里，例如MP2, CCSD, CAS等波函数含有多个行列式，其自然轨道没有严格的占据与空的概念，占据数一般也不是整数，也是[0,2]范围的小数。注意像MP2和CCSD的参考态RHF是单行列式的，但（请按/符号断句）/MP2里的一阶波函数/和/CCSD波函数/是多行列式的；而CASCI和CASSCF波函数本身就是多行列式的。
对一般的波函数而言，将电荷密度展开至原子基函数上（简便起见，省略上标NO）
写成矩阵形式即为
其中方阵的非对角元全是0，对角元。在此一般形式的波函数下，若自然轨道经过酉变换
中间的矩阵不再是对角矩阵，只是对称矩阵。因此对于一般形式的波函数，自然轨道是唯一的，正则轨道或其他类型的轨道不能充当自然轨道的角色。假设我们知道原子基密度矩阵（例如Gaussian的fchk文件里就有），如何求自然轨道占据数和自然轨道呢？
2.从密度矩阵求自然轨道
直接对角化矩阵是不行的，因为（1）自然轨道不是酉矩阵；（2）没法保证矩阵本征值的和等于总电子数。注意到我们有正交归一关系，我们可以给矩阵左右各乘一个
关于可阅读公众号本期另一文《的一些性质》。此矩阵的迹
便是总电子数，符合要求。可能有读者会有疑问，非得乘？事实上，仅考虑电子数的话，不止一种方式，如下通式均满足要求
但是，仅是对称矩阵。接着事情就很简单了，我们可以将这个对称矩阵对角化，
对比上述刚乘时的形式可以发现
则自然轨道系数矩阵为
在实际编程中求时需要舍弃接近零的值，即处理线性依赖。相应地，本征值得自己从大到小排序（MKL库函数输出是从小到大），取到自然分子轨道数目即止。若有本征值被舍弃，则的对应本征矢也应该舍弃，保证最后自然轨道系数矩阵的维度是基函数*自然轨道数。
注意，自然轨道数小于等于总轨道数。例如在CASCI和CASSCF方法中，若提供的密度矩阵是活性空间密度矩阵，则求出来的自然轨道数只能等于活性轨道数。若提供的密度矩阵是总密度矩阵，则自然轨道数等于总轨道数。
我们已经知道如何求自然轨道及其占据数，接着回到原有的情形：假设有一套普通的正交归一轨道（例如，Boys局域轨道，UNO轨道等等），它是自然轨道的酉变换，则对应的占据数矩阵为
可见，只需先求出变换关系，（可以调MKL库中解线性方程组的函数），然后做两次矩阵乘法即可得到。由于该轨道不是自然轨道，没有严格的占据数概念，其占据数矩阵不是对角的，在非对角元上也有值。我们可以将的对角元“近似地看成”第个轨道的“占据数”。假设的非对角元较小、对角元接近本征值，便可认为这套轨道与自然轨道较为接近，可以作为一种衡量接近自然轨道程度的指标。
最后，回到本文一开始的公式，假设我们现在将电荷密度展开在这组普通的正交归一轨道上
对比上文的，可以发现就是。
Acknowledgement
感谢清癯、zhigang、Acid、yuqiwang、暖云大师和食肉动物等人的审阅和建议。
Reference
Modern Quantum Chemistry, A. Szabo and N. S. Ostlund, p139. 《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》http://sobereva.com/403 http://gaussian.com/population
后记
本文使用mdnice的Chrome浏览器插件（支持Markdown和LaTex语法）写成。本篇为理论篇，在不确定的将来某天还有实战篇，先看看反响如何😂